function omega_hat = kernelEstimator(X, u, S, kernelType)
%KERNELESTIMATOR Estimates omega using Newey-West or QS kernel
%
%   omega_hat = kernelEstimator(X, u, S, kernelType)
%
%   Inputs:
%       X: T x k vector of regressor
%       u: T x 1 vector of residuals
%       S: bandwidth parameter
%       kernelType: 'NW' for Newey-West, 'QS' for Quadratic Spectral
%
%   Outputs:
%       omega_hat: estimated omega

z = X.*u;

T = size(z, 1);
k = size(z, 2);

omega_hat = zeros(k, k);

% Calculate Gamma_hat
for j = -(T-1):(T-1)
    Gamma_hat = zeros(k, k);

    v = j/S;
    if strcmp(kernelType, 'NW')
        ker = (1 - abs(v))*(abs(v) <= 1);
    elseif strcmp(kernelType, 'QS')
        x = 6*v/5;
        ker = 3*(sin(pi*x)/(pi*x) - cos(pi*x))/(pi*x)^2;
    else
        error('Invalid kernel type. Use "NW" or "QS".')
    end

    for t = max(1, j+1):min(T, T+j)
        Gamma_hat = Gamma_hat + z(t, :)'*z(t-j, :);
    end

    % Apply kernel function
    omega_hat = omega_hat + ker.*Gamma_hat;
end

    omega_hat = omega_hat./(T-k); % Small sample adjustment (Andrews, 1991).

    omega_hat = T.*((X'*X)\omega_hat/(X'*X));

end
